#### Performance legitimacy
#### v8: CPS R&R

library(lavaan)
library(lme4)
library(xtable)
library(texreg)
library(arm)
library(tidyr)
library(RColorBrewer)
library(rgdal)
library(psych)
library(mvtnorm)
library(plm)

# Setup and working directory
par(mfrow=c(1,1), mar=c(3,3,1,1), cex=0.9, las=1, mgp=c(1,0.4,0), tcl=-0.2)
WD = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(WD))
print( getwd() )

# Read data
sd = read.csv("replication_data.csv")
sd$Country = as.factor(sd$Country)
names(sd)

# Create democracy only datasets
sd$ClosedAut = ifelse(sd$Regime_VD==0, 1, 0)
sd$ClosedAut = ifelse(is.na(sd$Regime_VD), 1, sd$ClosedAut)
sd$Democ = ifelse(sd$Regime_VD>1, 1, 0)
sd$Democ = ifelse(is.na(sd$Regime_VD), 0, sd$Democ)

# Find first year of democracy
sd.d1 = sd[sd$Democ==1,]
first.yr.democ = data.frame(Country = levels(sd$Country), First_yr_democ = as.vector(
  by(sd.d1, sd.d1$Country, function(x) min(as.numeric(x$Year)))))
first.yr.democ[is.na(first.yr.democ$First_yr_democ), "First_yr_democ"] = 2020

# Edit first year for countries which later became closed autocracies
first.yr.democ[first.yr.democ$Country=="Thailand", "First_yr_democ"] = 2020
first.yr.democ[first.yr.democ$Country=="Slovakia", "First_yr_democ"] = 1995
first.yr.democ[first.yr.democ$Country=="Madagascar", "First_yr_democ"] = 2020
first.yr.democ[first.yr.democ$Country=="Libya", "First_yr_democ"] = 2020
first.yr.democ[first.yr.democ$Country=="Kyrgyzstan", "First_yr_democ"] = 2020
first.yr.democ[first.yr.democ$Country=="Kenya", "First_yr_democ"] = 2020
first.yr.democ[first.yr.democ$Country=="Togo", "First_yr_democ"] = 2020
sd.d.full = merge(sd, first.yr.democ, by="Country", all.x=TRUE)

# drop years before first year
sd.d.full = sd.d.full[sd.d.full$Year>=sd.d.full$First_yr_democ,] # dataset of democracies
supdem.d = sd.d.full[!is.na(sd.d.full$SupDem_z),]
supdem.d$Country = as.factor(as.character(supdem.d$Country))
# supdem.d = supdem.d[!(supdem.d$Country=="Bangladesh" & sd.d.full$Year==2008),] to remove Bangl2008 - closed autocracy


# Drop countries with less than 5 years' data
supdem.91 = supdem.d
drop.cnt = c("Ivory Coast", "Niger", "Suriname", "Trinidad and Tobago", "Tunisia", "Mauritius", "Sierra Leone", "Nigeria")
supdem.91 = supdem.91[!(supdem.91$Country %in% drop.cnt),]

# Create dataset of democracies using a stricter definition
supdem.80 = supdem.91
drop.cnt = c("Armenia", "Bangladesh", "Belarus", "Honduras", "Montenegro", "Mozambique", "Nicaragua", "Russia", "Sri Lanka", 
             "Ukraine", "Venezuela")
supdem.80 = supdem.80[!(supdem.80$Country %in% drop.cnt),]


## Full dataset of democracies

# Remove observations without 2nd lag of SupDem & 2nd lag of Satis
supdem.99.t1 = supdem.d[complete.cases(supdem.d[, c("Satis_m2", "SupDem_m2")]),]
supdem.99.t1$Country = as.factor(as.character(supdem.99.t1$Country))

# Create within and between measures of democracy and support
supdem.99.t1$SupDem_btw = ave(supdem.99.t1$SupDem_z, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$SupDem_within = supdem.99.t1$SupDem_z - supdem.99.t1$SupDem_btw
supdem.99.t1$SupDem_m1_btw = ave(supdem.99.t1$SupDem_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$SupDem_m1_within = supdem.99.t1$SupDem_m1 - supdem.99.t1$SupDem_m1_btw
supdem.99.t1$SupDem_m2_btw = ave(supdem.99.t1$SupDem_m2, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$SupDem_m2_within = supdem.99.t1$SupDem_m2 - supdem.99.t1$SupDem_m2_btw
supdem.99.t1$Satis_btw = ave(supdem.99.t1$Satis_z, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$Satis_within = supdem.99.t1$Satis_z - supdem.99.t1$Satis_btw
supdem.99.t1$Satis_m1_btw = ave(supdem.99.t1$Satis_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$Satis_m1_within = supdem.99.t1$Satis_m1 - supdem.99.t1$Satis_m1_btw
supdem.99.t1$Satis_m2_btw = ave(supdem.99.t1$Satis_m2, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$Satis_m2_within = supdem.99.t1$Satis_m2 - supdem.99.t1$Satis_m2_btw
supdem.99.t1$Poly_btw = ave(supdem.99.t1$Poly_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$Poly_within = supdem.99.t1$Poly_m1 - supdem.99.t1$Poly_btw
supdem.99.t1$Liberal_btw = ave(supdem.99.t1$Lib_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$Liberal_within = supdem.99.t1$Lib_m1 - supdem.99.t1$Liberal_btw
supdem.99.t1$GDP_btw = ave(supdem.99.t1$lnGDP_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$GDP_within = supdem.99.t1$lnGDP_m1 - supdem.99.t1$GDP_btw
supdem.99.t1$GDPgr_btw = ave(supdem.99.t1$GDP_grth_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$GDPgr_within = supdem.99.t1$GDP_grth_m1 - supdem.99.t1$GDPgr_btw
supdem.99.t1$Infl_btw = ave(supdem.99.t1$Inflat_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$Infl_within = supdem.99.t1$Inflat_m1 - supdem.99.t1$Infl_btw
supdem.99.t1$FreeViol_btw = ave(supdem.99.t1$FreeViol_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$FreeViol_within = supdem.99.t1$FreeViol_m1 - supdem.99.t1$FreeViol_btw
supdem.99.t1$ChgReg_btw = ave(supdem.99.t1$ChgReg_m1, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$ChgReg_within = supdem.99.t1$ChgReg_m1 - supdem.99.t1$ChgReg_btw
supdem.99.t1$ExecAppr_btw = ave(supdem.99.t1$ExecAppr_z, supdem.99.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t1$ExecAppr_within = supdem.99.t1$ExecAppr_z - supdem.99.t1$ExecAppr_btw

# Drop observations without lagged health, crime & corruption
supdem.99.t2 = supdem.99.t1[complete.cases(supdem.99.t1[, c("HAQI_m1", "ViolDeath_m1", "Corr_m1")]),]
supdem.99.t2$Country = as.factor(as.character(supdem.99.t2$Country))

# Create within and between measures of democracy and support
supdem.99.t2$SupDem_btw = ave(supdem.99.t2$SupDem_z, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$SupDem_within = supdem.99.t2$SupDem_z - supdem.99.t2$SupDem_btw
supdem.99.t2$SupDem_m1_btw = ave(supdem.99.t2$SupDem_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$SupDem_m1_within = supdem.99.t2$SupDem_m1 - supdem.99.t2$SupDem_m1_btw
supdem.99.t2$SupDem_m2_btw = ave(supdem.99.t2$SupDem_m2, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$SupDem_m2_within = supdem.99.t2$SupDem_m2 - supdem.99.t2$SupDem_m2_btw
supdem.99.t2$Satis_btw = ave(supdem.99.t2$Satis_z, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Satis_within = supdem.99.t2$Satis_z - supdem.99.t2$Satis_btw
supdem.99.t2$Satis_m1_btw = ave(supdem.99.t2$Satis_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Satis_m1_within = supdem.99.t2$Satis_m1 - supdem.99.t2$Satis_m1_btw
supdem.99.t2$Satis_m2_btw = ave(supdem.99.t2$Satis_m2, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Satis_m2_within = supdem.99.t2$Satis_m2 - supdem.99.t2$Satis_m2_btw
supdem.99.t2$Poly_btw = ave(supdem.99.t2$Poly_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Poly_within = supdem.99.t2$Poly_m1 - supdem.99.t2$Poly_btw
supdem.99.t2$Liberal_btw = ave(supdem.99.t2$Lib_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Liberal_within = supdem.99.t2$Lib_m1 - supdem.99.t2$Liberal_btw
supdem.99.t2$GDP_btw = ave(supdem.99.t2$lnGDP_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$GDP_within = supdem.99.t2$lnGDP_m1 - supdem.99.t2$GDP_btw
supdem.99.t2$GDPgr_btw = ave(supdem.99.t2$GDP_grth_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$GDPgr_within = supdem.99.t2$GDP_grth_m1 - supdem.99.t2$GDPgr_btw
supdem.99.t2$Infl_btw = ave(supdem.99.t2$Inflat_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Infl_within = supdem.99.t2$Inflat_m1 - supdem.99.t2$Infl_btw
supdem.99.t2$Unemp_btw = ave(supdem.99.t2$Unemp_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Unemp_within = supdem.99.t2$Unemp_m1 - supdem.99.t2$Unemp_btw
supdem.99.t2$Viol_btw = ave(supdem.99.t2$ViolDeath_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Viol_within = supdem.99.t2$ViolDeath_m1 - supdem.99.t2$Viol_btw
supdem.99.t2$InfMort_btw = ave(supdem.99.t2$InfMort_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$InfMort_within = supdem.99.t2$InfMort_m1 - supdem.99.t2$InfMort_btw
supdem.99.t2$Unemp_btw = ave(supdem.99.t2$Unemp_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Unemp_within = supdem.99.t2$Unemp_m1 - supdem.99.t2$Unemp_btw
supdem.99.t2$HAQI_btw = ave(supdem.99.t2$HAQI_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$HAQI_within = supdem.99.t2$HAQI_m1 - supdem.99.t2$HAQI_btw
supdem.99.t2$Corr_BCI_btw = ave(supdem.99.t2$Corr_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$Corr_BCI_within = supdem.99.t2$Corr_m1 - supdem.99.t2$Corr_BCI_btw
supdem.99.t2$FreeViol_btw = ave(supdem.99.t2$FreeViol_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$FreeViol_within = supdem.99.t2$FreeViol_m1 - supdem.99.t2$FreeViol_btw
supdem.99.t2$ChgReg_btw = ave(supdem.99.t2$ChgReg_m1, supdem.99.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.99.t2$ChgReg_within = supdem.99.t2$ChgReg_m1 - supdem.99.t2$ChgReg_btw

## 5+ years dataset of 91 democracies

# Drop observations without 2nd lag of SupDem & 2nd lag of Satis
supdem.91.t1 = supdem.91[complete.cases(supdem.91[, c("Satis_m2", "SupDem_m2")]),]
supdem.91.t1$Country = as.factor(as.character(supdem.91.t1$Country))

# Create within and between measures of democracy and support
supdem.91.t1$SupDem_btw = ave(supdem.91.t1$SupDem_z, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$SupDem_within = supdem.91.t1$SupDem_z - supdem.91.t1$SupDem_btw
supdem.91.t1$SupDem_m1_btw = ave(supdem.91.t1$SupDem_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$SupDem_m1_within = supdem.91.t1$SupDem_m1 - supdem.91.t1$SupDem_m1_btw
supdem.91.t1$SupDem_m2_btw = ave(supdem.91.t1$SupDem_m2, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$SupDem_m2_within = supdem.91.t1$SupDem_m2 - supdem.91.t1$SupDem_m2_btw
supdem.91.t1$Satis_btw = ave(supdem.91.t1$Satis_z, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$Satis_within = supdem.91.t1$Satis_z - supdem.91.t1$Satis_btw
supdem.91.t1$Satis_m1_btw = ave(supdem.91.t1$Satis_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$Satis_m1_within = supdem.91.t1$Satis_m1 - supdem.91.t1$Satis_m1_btw
supdem.91.t1$Satis_m2_btw = ave(supdem.91.t1$Satis_m2, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$Satis_m2_within = supdem.91.t1$Satis_m2 - supdem.91.t1$Satis_m2_btw
supdem.91.t1$Poly_btw = ave(supdem.91.t1$Poly_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$Poly_within = supdem.91.t1$Poly_m1 - supdem.91.t1$Poly_btw
supdem.91.t1$Liberal_btw = ave(supdem.91.t1$Lib_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$Liberal_within = supdem.91.t1$Lib_m1 - supdem.91.t1$Liberal_btw
supdem.91.t1$GDP_btw = ave(supdem.91.t1$lnGDP_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$GDP_within = supdem.91.t1$lnGDP_m1 - supdem.91.t1$GDP_btw
supdem.91.t1$GDPgr_btw = ave(supdem.91.t1$GDP_grth_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$GDPgr_within = supdem.91.t1$GDP_grth_m1 - supdem.91.t1$GDPgr_btw
supdem.91.t1$Infl_btw = ave(supdem.91.t1$Inflat_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$Infl_within = supdem.91.t1$Inflat_m1 - supdem.91.t1$Infl_btw
supdem.91.t1$FreeViol_btw = ave(supdem.91.t1$FreeViol_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$FreeViol_within = supdem.91.t1$FreeViol_m1 - supdem.91.t1$FreeViol_btw
supdem.91.t1$ChgReg_btw = ave(supdem.91.t1$ChgReg_m1, supdem.91.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t1$ChgReg_within = supdem.91.t1$ChgReg_m1 - supdem.91.t1$ChgReg_btw

# Drop observations without lagged health, crime & corruption
supdem.91.t2 = supdem.91.t1[complete.cases(supdem.91.t1[, c("HAQI_m1", "ViolDeath_m1", "Corr_m1")]),]
supdem.91.t2$Country = as.factor(as.character(supdem.91.t2$Country))

# Create within and between measures of democracy and support
supdem.91.t2$SupDem_btw = ave(supdem.91.t2$SupDem_z, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$SupDem_within = supdem.91.t2$SupDem_z - supdem.91.t2$SupDem_btw
supdem.91.t2$SupDem_m1_btw = ave(supdem.91.t2$SupDem_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$SupDem_m1_within = supdem.91.t2$SupDem_m1 - supdem.91.t2$SupDem_m1_btw
supdem.91.t2$SupDem_m2_btw = ave(supdem.91.t2$SupDem_m2, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$SupDem_m2_within = supdem.91.t2$SupDem_m2 - supdem.91.t2$SupDem_m2_btw
supdem.91.t2$Satis_btw = ave(supdem.91.t2$Satis_z, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Satis_within = supdem.91.t2$Satis_z - supdem.91.t2$Satis_btw
supdem.91.t2$Satis_m1_btw = ave(supdem.91.t2$Satis_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Satis_m1_within = supdem.91.t2$Satis_m1 - supdem.91.t2$Satis_m1_btw
supdem.91.t2$Satis_m2_btw = ave(supdem.91.t2$Satis_m2, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Satis_m2_within = supdem.91.t2$Satis_m2 - supdem.91.t2$Satis_m2_btw
supdem.91.t2$Poly_btw = ave(supdem.91.t2$Poly_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Poly_within = supdem.91.t2$Poly_m1 - supdem.91.t2$Poly_btw
supdem.91.t2$Liberal_btw = ave(supdem.91.t2$Lib_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Liberal_within = supdem.91.t2$Lib_m1 - supdem.91.t2$Liberal_btw
supdem.91.t2$GDP_btw = ave(supdem.91.t2$lnGDP_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$GDP_within = supdem.91.t2$lnGDP_m1 - supdem.91.t2$GDP_btw
supdem.91.t2$GDPgr_btw = ave(supdem.91.t2$GDP_grth_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$GDPgr_within = supdem.91.t2$GDP_grth_m1 - supdem.91.t2$GDPgr_btw
supdem.91.t2$Infl_btw = ave(supdem.91.t2$Inflat_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Infl_within = supdem.91.t2$Inflat_m1 - supdem.91.t2$Infl_btw
supdem.91.t2$Unemp_btw = ave(supdem.91.t2$Unemp_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Unemp_within = supdem.91.t2$Unemp_m1 - supdem.91.t2$Unemp_btw
supdem.91.t2$Viol_btw = ave(supdem.91.t2$ViolDeath_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Viol_within = supdem.91.t2$ViolDeath_m1 - supdem.91.t2$Viol_btw
supdem.91.t2$InfMort_btw = ave(supdem.91.t2$InfMort_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$InfMort_within = supdem.91.t2$InfMort_m1 - supdem.91.t2$InfMort_btw
supdem.91.t2$Unemp_btw = ave(supdem.91.t2$Unemp_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Unemp_within = supdem.91.t2$Unemp_m1 - supdem.91.t2$Unemp_btw
supdem.91.t2$HAQI_btw = ave(supdem.91.t2$HAQI_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$HAQI_within = supdem.91.t2$HAQI_m1 - supdem.91.t2$HAQI_btw
supdem.91.t2$Corr_BCI_btw = ave(supdem.91.t2$Corr_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$Corr_BCI_within = supdem.91.t2$Corr_m1 - supdem.91.t2$Corr_BCI_btw
supdem.91.t2$FreeViol_btw = ave(supdem.91.t2$FreeViol_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$FreeViol_within = supdem.91.t2$FreeViol_m1 - supdem.91.t2$FreeViol_btw
supdem.91.t2$ChgReg_btw = ave(supdem.91.t2$ChgReg_m1, supdem.91.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.91.t2$ChgReg_within = supdem.91.t2$ChgReg_m1 - supdem.91.t2$ChgReg_btw

## Dataset of 80 stringent democracies

# Drop observations without 2nd lag of SupDem & 2nd lag of Satis
supdem.80.t1 = supdem.80[complete.cases(supdem.80[, c("Satis_m2", "SupDem_m2")]),]
supdem.80.t1$Country = as.factor(as.character(supdem.80.t1$Country))

# Create within and between measures of democracy and support
supdem.80.t1$SupDem_btw = ave(supdem.80.t1$SupDem_z, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$SupDem_within = supdem.80.t1$SupDem_z - supdem.80.t1$SupDem_btw
supdem.80.t1$SupDem_m1_btw = ave(supdem.80.t1$SupDem_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$SupDem_m1_within = supdem.80.t1$SupDem_m1 - supdem.80.t1$SupDem_m1_btw
supdem.80.t1$SupDem_m2_btw = ave(supdem.80.t1$SupDem_m2, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$SupDem_m2_within = supdem.80.t1$SupDem_m2 - supdem.80.t1$SupDem_m2_btw
supdem.80.t1$Satis_btw = ave(supdem.80.t1$Satis_z, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$Satis_within = supdem.80.t1$Satis_z - supdem.80.t1$Satis_btw
supdem.80.t1$Satis_m1_btw = ave(supdem.80.t1$Satis_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$Satis_m1_within = supdem.80.t1$Satis_m1 - supdem.80.t1$Satis_m1_btw
supdem.80.t1$Satis_m2_btw = ave(supdem.80.t1$Satis_m2, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$Satis_m2_within = supdem.80.t1$Satis_m2 - supdem.80.t1$Satis_m2_btw
supdem.80.t1$Poly_btw = ave(supdem.80.t1$Poly_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$Poly_within = supdem.80.t1$Poly_m1 - supdem.80.t1$Poly_btw
supdem.80.t1$Liberal_btw = ave(supdem.80.t1$Lib_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$Liberal_within = supdem.80.t1$Lib_m1 - supdem.80.t1$Liberal_btw
supdem.80.t1$GDP_btw = ave(supdem.80.t1$lnGDP_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$GDP_within = supdem.80.t1$lnGDP_m1 - supdem.80.t1$GDP_btw
supdem.80.t1$GDPgr_btw = ave(supdem.80.t1$GDP_grth_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$GDPgr_within = supdem.80.t1$GDP_grth_m1 - supdem.80.t1$GDPgr_btw
supdem.80.t1$Infl_btw = ave(supdem.80.t1$Inflat_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$Infl_within = supdem.80.t1$Inflat_m1 - supdem.80.t1$Infl_btw
supdem.80.t1$Viol_btw = ave(supdem.80.t1$ViolDeath_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$Viol_within = supdem.80.t1$ViolDeath_m1 - supdem.80.t1$Viol_btw
supdem.80.t1$InfMort_btw = ave(supdem.80.t1$InfMort_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$InfMort_within = supdem.80.t1$InfMort_m1 - supdem.80.t1$InfMort_btw
supdem.80.t1$Unemp_btw = ave(supdem.80.t1$Unemp_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t1$Unemp_within = supdem.80.t1$Unemp_m1 - supdem.80.t1$Unemp_btw
supdem.80.t1$HAQI_btw = ave(supdem.80.t1$HAQI_m1, supdem.80.t1$Country, FUN=function(x) mean(x, na.rm=TRUE))
supdem.80.t1$HAQI_within = supdem.80.t1$HAQI_m1 - supdem.80.t1$HAQI_btw

# Drop observations without lagged health, crime & corruption
supdem.80.t2 = supdem.80.t1[complete.cases(supdem.80.t1[, c("HAQI_m1", "ViolDeath_m1", "Corr_m1")]),]
supdem.80.t2$Country = as.factor(as.character(supdem.80.t2$Country))

# Create within and between measures of democracy and support
supdem.80.t2$SupDem_btw = ave(supdem.80.t2$SupDem_z, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$SupDem_within = supdem.80.t2$SupDem_z - supdem.80.t2$SupDem_btw
supdem.80.t2$SupDem_m1_btw = ave(supdem.80.t2$SupDem_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$SupDem_m1_within = supdem.80.t2$SupDem_m1 - supdem.80.t2$SupDem_m1_btw
supdem.80.t2$SupDem_m2_btw = ave(supdem.80.t2$SupDem_m2, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$SupDem_m2_within = supdem.80.t2$SupDem_m2 - supdem.80.t2$SupDem_m2_btw
supdem.80.t2$Satis_btw = ave(supdem.80.t2$Satis_z, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Satis_within = supdem.80.t2$Satis_z - supdem.80.t2$Satis_btw
supdem.80.t2$Satis_m1_btw = ave(supdem.80.t2$Satis_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Satis_m1_within = supdem.80.t2$Satis_m1 - supdem.80.t2$Satis_m1_btw
supdem.80.t2$Satis_m2_btw = ave(supdem.80.t2$Satis_m2, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Satis_m2_within = supdem.80.t2$Satis_m2 - supdem.80.t2$Satis_m2_btw
supdem.80.t2$Poly_btw = ave(supdem.80.t2$Poly_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Poly_within = supdem.80.t2$Poly_m1 - supdem.80.t2$Poly_btw
supdem.80.t2$Liberal_btw = ave(supdem.80.t2$Lib_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Liberal_within = supdem.80.t2$Lib_m1 - supdem.80.t2$Liberal_btw
supdem.80.t2$GDP_btw = ave(supdem.80.t2$lnGDP_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$GDP_within = supdem.80.t2$lnGDP_m1 - supdem.80.t2$GDP_btw
supdem.80.t2$GDPgr_btw = ave(supdem.80.t2$GDP_grth_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$GDPgr_within = supdem.80.t2$GDP_grth_m1 - supdem.80.t2$GDPgr_btw
supdem.80.t2$Infl_btw = ave(supdem.80.t2$Inflat_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Infl_within = supdem.80.t2$Inflat_m1 - supdem.80.t2$Infl_btw
supdem.80.t2$Unemp_btw = ave(supdem.80.t2$Unemp_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Unemp_within = supdem.80.t2$Unemp_m1 - supdem.80.t2$Unemp_btw
supdem.80.t2$Viol_btw = ave(supdem.80.t2$ViolDeath_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Viol_within = supdem.80.t2$ViolDeath_m1 - supdem.80.t2$Viol_btw
supdem.80.t2$InfMort_btw = ave(supdem.80.t2$InfMort_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$InfMort_within = supdem.80.t2$InfMort_m1 - supdem.80.t2$InfMort_btw
supdem.80.t2$Unemp_btw = ave(supdem.80.t2$Unemp_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Unemp_within = supdem.80.t2$Unemp_m1 - supdem.80.t2$Unemp_btw
supdem.80.t2$HAQI_btw = ave(supdem.80.t2$HAQI_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$HAQI_within = supdem.80.t2$HAQI_m1 - supdem.80.t2$HAQI_btw
supdem.80.t2$Corr_BCI_btw = ave(supdem.80.t2$Corr_m1, supdem.80.t2$Country, FUN=function(x) mean(x, na.rm=FALSE))
supdem.80.t2$Corr_BCI_within = supdem.80.t2$Corr_m1 - supdem.80.t2$Corr_BCI_btw


## Create panel datasets

# panel dataset of democracies, with countries with less than 3 years support data dropped.
sd.plm.d = pdata.frame(supdem.d, index = c("Country", "Year"))
table(index(sd.plm.d), useNA = "ifany")

# panel dataset of democracies with 5+ years data, with support, satis, murder, ihme health, corrup data
sd.plm.91.t2 = pdata.frame(supdem.91.t2, index = c("Country", "Year"))



### Descriptive figures

## Dem support and satisfaction

col1 = rgb(1, 0.55, 0, 1)
col2 = rgb(0, 0.45, 0.7, 1)
cnts = levels(supdem.91.t1$Country)
n.cnts = length(cnts)

pdf("fig_1.pdf", width=6, height=7.5)
par(mfrow=c(8,6), mar=c(1,1.2,0.2,0.2), tcl=-0.15, cex=0.75)
for(i in 1:48) {
  plot(x=supdem.91.t1[supdem.91.t1$Country==cnts[i], "Year"],
       y=supdem.91.t1[supdem.91.t1$Country==cnts[i], "SupDem_z"],
       type="n", xlim=c(1988, 2018), ylim=c(-2.6, 3.5), xlab="", ylab="", main="", axes=FALSE)
  abline(h=c(-2,0,2), lty=3, lwd=0.5, col=grey(0.5))
  abline(v=c(1990, 2000, 2010), lty=3, lwd=0.5, col=grey(0.5))
  lines(x=supdem.91.t1[supdem.91.t1$Country==cnts[i], "Year"],
        y=supdem.91.t1[supdem.91.t1$Country==cnts[i], "SupDem_z"],
        col=col1, lwd=1.5)
  lines(x=supdem.91.t1[supdem.91.t1$Country==cnts[i], "Year"],
        y=supdem.91.t1[supdem.91.t1$Country==cnts[i], "Satis_z"],
        col=col2, lwd=1.5)
  text(cnts[i], x=1988, y=3.2, cex=0.7, adj=0)
  axis(side=1, at=c(1990,2000,2010), labels=c(1990,2000,2010), mgp=c(0,-0.1,0), cex.axis=0.7)
  axis(side=2, at=c(-2,0,2), labels=c(-2,0,2), mgp=c(0,0.1,0), cex.axis=0.7)
  box()
}
text("Support", x=1996, y=supdem.91.t1[supdem.91.t1$Country==cnts[48] & supdem.91.t1$Year==2008, "SupDem_z"] - 0.5, 
     cex=0.65, adj=0, col=col1)
text("Satisfaction", x=1997, y=supdem.91.t1[supdem.91.t1$Country==cnts[48] & supdem.91.t1$Year==2012, "Satis_z"] + 0.3, 
     cex=0.65, adj=0, col=col2)
dev.off()


pdf("fig_2.pdf", width=6, height=7.5)
par(mfrow=c(8,6), mar=c(1,1.2,0.2,0.2), tcl=-0.15, cex=0.75)
for(i in 49:n.cnts) {
  plot(x=supdem.91.t1[supdem.91.t1$Country==cnts[i], "Year"],
       y=supdem.91.t1[supdem.91.t1$Country==cnts[i], "SupDem_z"],
       type="n", xlim=c(1988, 2018), ylim=c(-2.6, 3.5), xlab="", ylab="", main="", axes=FALSE)
  abline(h=c(-2,0,2), lty=3, lwd=0.5, col=grey(0.5))
  abline(v=c(1990, 2000, 2010), lty=3, lwd=0.5, col=grey(0.5))
  lines(x=supdem.91.t1[supdem.91.t1$Country==cnts[i], "Year"],
        y=supdem.91.t1[supdem.91.t1$Country==cnts[i], "SupDem_z"],
        col=col1, lwd=1.5)
  lines(x=supdem.91.t1[supdem.91.t1$Country==cnts[i], "Year"],
        y=supdem.91.t1[supdem.91.t1$Country==cnts[i], "Satis_z"],
        col=col2, lwd=1.5)
  text(cnts[i], x=1988, y=3.2, cex=0.7, adj=0)
  axis(side=1, at=c(1990,2000,2010), labels=c(1990,2000,2010), mgp=c(0,-0.1,0), cex.axis=0.7)
  axis(side=2, at=c(-2,0,2), labels=c(-2,0,2), mgp=c(0,0.1,0), cex.axis=0.7)
  box()
}
text("Support", x=1992, y=supdem.91.t1[supdem.91.t1$Country==cnts[length(cnts)] & supdem.91.t1$Year==2005, "SupDem_z"] + 0.3, 
     cex=0.65, adj=0, col=col1)
text("Satisfaction", x=1993, y=supdem.91.t1[supdem.91.t1$Country==cnts[length(cnts)] & supdem.91.t1$Year==2005, "Satis_z"] - 0.6, 
     cex=0.65, adj=0, col=col2)
dev.off()


## Distributions of countries across regime types

cnts = levels(supdem.91.t1$Country)
supdem.mat = spread(supdem.91.t1[,c(1,2,11)], key=Year, value=Regime_VD)
regime.mat = t(as.matrix(supdem.mat[,2:30]))

reg.pal = c("#ca0020", "#f4a582", "#92c5de", "#0571b0")
n.cnt = length(levels(supdem.91.t1$Country))
pdf("fig_s1.pdf", height=10, width=6)
par(mfrow=c(1,1), mar=c(2.5, 9.3, 1.5, 0.5), cex=0.6, tcl=-0.2, mgp=c(1,0.2,0), las=1)
image(x=1990:2018, y=1:n.cnt, z=regime.mat[,n.cnt:1], col=reg.pal, yaxt="n", xaxt="n", xlab="Year",
      ylab="", cex.lab=1.2, useRaster=FALSE)
abline(v=1990:2018-0.5, col=rgb(0,0,0,0.5), lty=1, lwd=0.5)
abline(h=(1:(n.cnt+1))-0.5, col=rgb(0,0,0,0.5), lty=1, lwd=0.5)
axis(1, at=(c(1990,1995,2000,2005,2010,2015))-0.5, labels=c(1990,1995,2000,2005,2010,2015), 
     tick=TRUE, cex.axis=0.78, mgp=c(1,0.1,0))
axis(2, at=(1:n.cnt)+0.1, tick=FALSE, labels=cnts[n.cnt:1], cex.axis=1)
image(x=1990:2018, y=1:n.cnt, z=regime.mat[,n.cnt:1], col=reg.pal, add=TRUE, useRaster=FALSE)
abline(v=1990:2018-0.5, col=rgb(1,1,1,0.5), lty=1, lwd=0.5)
abline(h=1:(n.cnt+1)-0.5, col=rgb(1,1,1,0.5), lty=1, lwd=0.5)
box()
dev.off()



### Models


## Correlations between various opinion series

summary(opin.test.1 <- plm(SupDem_z ~ Satis_z, effect="indiv", model="within", sd.plm.d))
summary(opin.test.2 <- plm(SupDem_z ~ ExecAppr_z, effect="indiv", model="within", sd.plm.d))
summary(opin.test.3 <- plm(SupDem_z ~ ExecAppr_z + Satis_z, effect="indiv", model="within", sd.plm.d))

opin.mods = list(opin.test.1, opin.test.2, opin.test.3)
sapply(opin.mods, function(x) round(sigma(x), 3)) # residual standard error
vcm = list()
vcm = lapply(opin.mods, function(x) plm::vcovHC(x, method="white1", cluster="group"))
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(opin.mods)) { pval[[i]] = pt(abs(coef(opin.mods[[i]]) / rse[[i]]), df=1000, lower.tail=FALSE) * 2 }
screenreg(opin.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE)

varnames = c("Democratic satisfaction", "Executive approval")

texreg(opin.mods, file="tab_s8.txt", override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, 
       leading.zero=FALSE, stars=0.05, custom.coef.names=varnames, booktabs=TRUE, dcolumn=TRUE,
       caption="Table S8: Fixed effects models of democratic support, satisfaction, and executive approval")


## Structural equation models, 91-democracy sample

# Basic model
sem.mod.eco <- 
'
SupDem_within ~ SupDem_m1_within + SupDem_m2_within + Satis_within + Poly_within + Liberal_within 
                + GDP_within + GDPgr_within + Infl_within
Satis_within ~ Satis_m1_within + Satis_m2_within + SupDem_within + Poly_within + Liberal_within 
                + GDP_within + GDPgr_within + Infl_within
'
sem.fit.eco.91 = sem(model=sem.mod.eco, data=supdem.91.t1, estimator="MLR")
summary(sem.fit.eco.91, nd=3)


# Full model
sem.mod.full <- 
  '
SupDem_within ~ SupDem_m1_within + SupDem_m2_within + Satis_within + Poly_within + Liberal_within 
                + GDP_within + GDPgr_within + Infl_within + HAQI_within + Viol_within + Corr_BCI_within
Satis_within ~ Satis_m1_within + Satis_m2_within + SupDem_within + Poly_within + Liberal_within 
                + GDP_within + GDPgr_within + Infl_within + HAQI_within + Viol_within + Corr_BCI_within
'
sem.fit.full.91 = sem(model=sem.mod.full, data=supdem.91.t2, estimator="MLR")
summary(sem.fit.full.91, nd=3)


# Extra covariates -- infant mortality and unemployment (for appendix)
sem.mod.addcov <- 
  '
SupDem_within ~ SupDem_m1_within + SupDem_m2_within + Satis_within + Poly_within + Liberal_within 
                + GDP_within + GDPgr_within + Infl_within + HAQI_within + Viol_within + Corr_BCI_within
                + InfMort_within + Unemp_within
Satis_within ~ Satis_m1_within + Satis_m2_within + SupDem_within + Poly_within + Liberal_within 
                + GDP_within + GDPgr_within + Infl_within + HAQI_within + Viol_within + Corr_BCI_within
                + InfMort_within + Unemp_within
'
sem.fit.addcov.91 = sem(model=sem.mod.addcov, data=supdem.91.t2, estimator="MLR")
summary(sem.fit.addcov.91, nd=3)


# Extra covariates -- regime change and polit violence (for appendix)
sem.mod.addcov2 <- 
  '
SupDem_within ~ SupDem_m1_within + SupDem_m2_within + Satis_within + Poly_within + Liberal_within 
                + GDP_within + GDPgr_within + Infl_within + HAQI_within + Viol_within + Corr_BCI_within
                + FreeViol_within + ChgReg_within
Satis_within ~ Satis_m1_within + Satis_m2_within + SupDem_within + Poly_within + Liberal_within 
                + GDP_within + GDPgr_within + Infl_within + HAQI_within + Viol_within + Corr_BCI_within
                + FreeViol_within + ChgReg_within
'
sem.fit.addcov2.91 = sem(model=sem.mod.addcov2, data=supdem.91.t2, estimator="MLR")
summary(sem.fit.addcov2.91, nd=3)


## Structural equation models with stringent 80-democracy sample

# Basic model
sem.fit.eco.80 = sem(model=sem.mod.eco, data=supdem.80.t1, estimator="MLR")
summary(sem.fit.eco.80, nd=3)

# Full model
sem.fit.full.80 = sem(model=sem.mod.full, data=supdem.80.t2, estimator="MLR")
summary(sem.fit.full.80, nd=3)


## Structural equation models with full 99-state sample

# Basic model
sem.fit.eco.99 = sem(model=sem.mod.eco, data=supdem.99.t1, estimator="MLR")
summary(sem.fit.eco.99, nd=3)

# Full model
sem.fit.full.99 = sem(model=sem.mod.full, data=supdem.99.t2, estimator="MLR")
summary(sem.fit.full.99, nd=3)


## FE models without cross-effects

sup.mod.2 = lm(SupDem_within ~ 0 + SupDem_m1_within + SupDem_m2_within + Poly_within + Liberal_within + GDP_within 
               + GDPgr_within + Infl_within + HAQI_within + Viol_within + Corr_BCI_within, supdem.91.t2)
sat.mod.2 = lm(Satis_within ~ 0 + Satis_m1_within + Satis_m2_within + Poly_within + Liberal_within + GDP_within 
               + GDPgr_within + Infl_within + HAQI_within + Viol_within + Corr_BCI_within, supdem.91.t2)
screenreg(list(sat.mod.2, sup.mod.2), digits=3)

fe1.mods = list(sup.mod.2, sat.mod.2)
screenreg(fe1.mods, digits=3)
sapply(fe1.mods, function(x) round(sigma(x), 3)) # residual standard error
vcm = list()
vcm = lapply(fe1.mods, function(x) plm::vcovHC(x, method="white1", cluster="group"))
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(fe1.mods)) { pval[[i]] = pt(abs(coef(fe1.mods[[i]]) / rse[[i]]), df=1600, lower.tail=FALSE) * 2 }
screenreg(fe1.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, leading.zero=FALSE)

varnames = c("Democratic support t-1", "Democratic support t-2", "Electoral democracy t-1", 
             "Liberalism t-1", "Log GDP/capita t-1", "GDP growth t-1", "Inflation rate t-1", 
             "Healthcare access & quality t-1", "Log rate of interpersonal violence t-1", 
             "Corruption index", "Democratic satisfaction t-1", "Democratic satisfaction t-2")
modnames = c("Sup1", "Sat1")

texreg(fe1.mods, file="tab_s7.txt", override.se=rse, override.pvalues=pval, digits=3, single.row=FALSE, 
       leading.zero=FALSE, stars=0.05, custom.coef.names=varnames, custom.model.names=modnames, booktabs=TRUE, dcolumn=TRUE,
       caption="Table S7: Separate fixed effects models of democratic support and satisfaction")


## Create SEM table 

capture.output(summary(sem.fit.eco.91, nd=3), file="tab_1.1.txt")
capture.output(summary(sem.fit.full.91, nd=3), file="tab_1.2.txt")

capture.output(summary(sem.fit.eco.80, nd=3), file="tab_s4.1.txt")
capture.output(summary(sem.fit.full.80, nd=3), file="tab_s4.2.txt")

capture.output(summary(sem.fit.eco.99, nd=3), file="tab_s3.1.txt")
capture.output(summary(sem.fit.full.99, nd=3), file="tab_s3.2.txt")

capture.output(summary(sem.fit.addcov.91, nd=3), file="tab_s5.txt")
capture.output(summary(sem.fit.addcov2.91, nd=3), file="tab_s6.txt")


#### Other tables

## Table of descriptives

descrip.var = c("SupDem_within", "Satis_within", "Poly_within", "Liberal_within", "GDP_within", "GDPgr_within", "Infl_within",
                "HAQI_within", "Viol_within", "Corr_BCI_within")
descrip.tab = matrix(NA, nrow=length(descrip.var), ncol=4)
descrip.tab[,1] = round(sapply(supdem.91.t2[,descrip.var], mean), 2)
descrip.tab[,2] = round(sapply(supdem.91.t2[,descrip.var], sd), 2)
descrip.tab[,3] = round(sapply(supdem.91.t2[,descrip.var], min), 2)
descrip.tab[,4] = round(sapply(supdem.91.t2[,descrip.var], max), 2)
colnames(descrip.tab) = c("Mean", "SD", "Min", "Max")
row.names(descrip.tab) = descrip.var
print(xtable(descrip.tab, digits=2), file="tab_s1.txt")


## Tests of number of lags

# support
sup.adl1.0 = lm(SupDem_z ~ Satis_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1 + Country, supdem.91.t2)
sup.adl1.1 = lm(SupDem_z ~ SupDem_m1 + Satis_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1 + Country, supdem.91.t2)
sup.adl1.2 = lm(SupDem_z ~ SupDem_m1 + SupDem_m2 + Satis_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1 + Country, supdem.91.t2)
sup.mods = list(sup.adl1.0, sup.adl1.1, sup.adl1.2)
screenreg(sup.mods, single.row=TRUE, digits=3)
sup.lag.tab = sapply(sup.mods, extractAIC) # 2 lags optimal
dimnames(sup.lag.tab) = list(c("df", "AIC"), c("0lag", "1lag", "2lag"))
round(sup.lag.tab, 1)

# satis
sat.adl1.0 = lm(Satis_z ~ SupDem_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1 + Country, supdem.91.t2)
sat.adl1.1 = lm(Satis_z ~ Satis_m1 + SupDem_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1 + Country, supdem.91.t2)
sat.adl1.2 = lm(Satis_z ~ Satis_m1 + Satis_m2 + SupDem_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1 + Country, supdem.91.t2)
sat.mods = list(sat.adl1.0, sat.adl1.1, sat.adl1.2)
screenreg(sat.mods, single.row=TRUE, digits=3)
sat.lag.tab = sapply(sat.mods, extractAIC) # 2 lags optimal
dimnames(sat.lag.tab) = list(c("df", "AIC"), c("0lag", "1lag", "2lag"))
round(sat.lag.tab, 1)

lag.tab = rbind(sup.lag.tab, sat.lag.tab)
print(xtable(lag.tab, digits=1), file="tab_s2a.txt")


## Autocorrelation tests

sup.adl1.0 = plm(SupDem_z ~ Satis_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1, sd.plm.91.t2, model="within")
sup.adl1.1 = plm(SupDem_z ~ SupDem_m1 + Satis_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1, sd.plm.91.t2, model="within")
sup.adl1.2 = plm(SupDem_z ~ SupDem_m1 + SupDem_m2 + Satis_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1, sd.plm.91.t2, model="within")
sup.mods = list(sup.adl1.0, sup.adl1.1, sup.adl1.2)
sapply(sup.mods, function(x) round(pwartest(x)[["p.value"]], 3)) # Wooldridge test p-val

sat.adl1.0 = plm(Satis_z ~ SupDem_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1, sd.plm.91.t2, model="within")
sat.adl1.1 = plm(Satis_z ~ Satis_m1 + SupDem_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1, sd.plm.91.t2, model="within")
sat.adl1.2 = plm(Satis_z ~ Satis_m1 + Satis_m2 + SupDem_z + Poly_m1 + Lib_m1 + lnGDP_m1 + GDP_grth_m1 + Inflat_m1
                + HAQI_m1 + ViolDeath_m1 + Corr_m1, sd.plm.91.t2, model="within")
sat.mods = list(sat.adl1.0, sat.adl1.1, sat.adl1.2)
sapply(sat.mods, function(x) round(pwartest(x)[["p.value"]], 3)) # Wooldridge test p-val




### Simulated Effects Plots 


## GDP growth

mod = sem.fit.full.91
round(coef(sem.fit.full.91), 2)
n.coef.sup =  length(coef(mod)) / 2 - 1
n.coef.sat =  length(coef(mod)) / 2 - 1
n.coef = length(coef(mod)) - 2
n.yrs = 121
n.burn = 100

summary(supdem.91.t1$GDPgr_within)
sd(supdem.91.t1$GDPgr_within, na.rm=TRUE)
chg.grth = 0.05

# support simulation

pred.data.sup = data.frame(DemSup_m1=0, DemSup_m2=0, DemSat=0, 
                           Poly=0, Lib=0, GDP=0, 
                           GDPgr=c(rep(0, n.burn), rep(chg.grth, n.yrs-n.burn)), 
                           Infl=0, HAQI=0, Viol=0, Corr=0,
                           DemSat_m1=0, DemSat_m2=0, DemSup=0)

sup.mod.nam = c("DemSup_m1", "DemSup_m2", "DemSat", "Poly", "Lib", "GDP", "GDPgr", "Infl", "HAQI", "Viol", "Corr")
sat.mod.nam = c("DemSat_m1", "DemSat_m2", "DemSup", "Poly", "Lib", "GDP", "GDPgr", "Infl", "HAQI", "Viol", "Corr")

n.sims = 1e3
sim.data = array(as.matrix(pred.data.sup), dim=c(dim(pred.data.sup), n.sims))
mod.sims = mvtnorm::rmvnorm(n=n.sims, mean=coef(mod), sigma=vcov(mod))

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    sim.data[i, match("DemSat", names(pred.data.sup)), k] =
      mod.sims[k,(n.coef.sat+1):n.coef] %*% as.matrix(sim.data[i, match(sat.mod.nam, names(pred.data.sup)), k])
    sim.data[i+1, match("DemSat_m1", names(pred.data.sup)), k] = sim.data[i, match("DemSat", names(pred.data.sup)), k]
    sim.data[i+1, match("DemSat_m2", names(pred.data.sup)), k] = sim.data[i, match("DemSat_m1", names(pred.data.sup)), k]
    sim.data[i, match("DemSup", names(pred.data.sup)), k] = 
      mod.sims[k,1:n.coef.sup] %*% as.matrix(sim.data[i, match(sup.mod.nam, names(pred.data.sup)), k]) 
    sim.data[i+1, match("DemSup_m1", names(pred.data.sup)), k] = sim.data[i, match("DemSup", names(pred.data.sup)), k]
    sim.data[i+1, match("DemSup_m2", names(pred.data.sup)), k] = sim.data[i, match("DemSup_m1", names(pred.data.sup)), k]
  }
}

pe.sup = apply(sim.data[, match("DemSup", names(pred.data.sup)), ], 1, mean, na.rm=TRUE)
pe.sd = apply(sim.data[, match("DemSup", names(pred.data.sup)), ], 1, sd, na.rm=TRUE)
u95.sup = pe.sup + 1.96*pe.sd
l95.sup = pe.sup - 1.96*pe.sd

yr.st = n.burn - 1

# satisfaction simulation

pred.data.sat = data.frame(DemSup_m1=0, DemSup_m2=0, DemSat=0, 
                           Poly=0, Lib=0, GDP=0, 
                           GDPgr=c(rep(0, n.burn), rep(chg.grth, n.yrs-n.burn)), 
                           Infl=0, HAQI=0, Viol=0, Corr=0,
                           DemSat_m1=0, DemSat_m2=0, DemSup=0)

n.sims = 1e3
sim.data = array(as.matrix(pred.data.sat), dim=c(dim(pred.data.sat), n.sims))

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    sim.data[i, match("DemSup", names(pred.data.sat)), k] = 
      mod.sims[k, 1:n.coef.sup] %*% as.matrix(sim.data[i, match(sup.mod.nam, names(pred.data.sat)), k]) 
    sim.data[i+1, match("DemSup_m1", names(pred.data.sat)), k] = sim.data[i, match("DemSup", names(pred.data.sat)), k]
    sim.data[i+1, match("DemSup_m2", names(pred.data.sat)), k] = sim.data[i, match("DemSup_m1", names(pred.data.sat)), k]
    sim.data[i, match("DemSat", names(pred.data.sat)), k] = 
      mod.sims[k,(n.coef.sat+1):n.coef] %*% as.matrix(sim.data[i, match(sat.mod.nam, names(pred.data.sat)), k]) 
    sim.data[i+1, match("DemSat_m1", names(pred.data.sat)), k] = sim.data[i, match("DemSat", names(pred.data.sat)), k]
    sim.data[i+1, match("DemSat_m2", names(pred.data.sat)), k] = sim.data[i, match("DemSat_m1", names(pred.data.sat)), k]
  }
}

pe.sat = apply(sim.data[, match("DemSat", names(pred.data.sat)), ], 1, mean, na.rm=TRUE)
pe.sd = apply(sim.data[, match("DemSat", names(pred.data.sat)), ], 1, sd, na.rm=TRUE)
u95.sat = pe.sat + 1.96*pe.sd
l95.sat = pe.sat - 1.96*pe.sd

# plot

pdf("fig_3a.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe.sup[yr.st:(n.yrs-1)], type="n", ylim=c(-0.2, 0.2), 
     xaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(110, 120, by=10), lty=3, lwd=0.75, col=grey(0.5))
abline(h=c(-0.2, -0.1, 0, 0.1, 0.2), lty=3, lwd=0.75, col=grey(0.5))
abline(v=100, lty=2, lwd=0.75, col=grey(0.5))
axis(side=1, at=seq(100, 120, by=10), labels=seq(0, 20, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-0.2, -0.1, 0, 0.1, 0.2), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe.sup[yr.st:(n.yrs-1)], col=rgb(0,.4,.5,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(0,.4,.5,.5), border=NA, 
        y=c(u95.sup[yr.st:(n.yrs-1)], l95.sup[(n.yrs-1):yr.st]))
lines(x=yr.st:(n.yrs-1), y=pe.sat[yr.st:(n.yrs-1)], col=rgb(1,.55,0,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(1,.55,0,.5), border=NA, 
        y=c(u95.sat[yr.st:(n.yrs-1)], l95.sat[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Democratic support / satisfaction", cex=0.9)
text(x=110, y=pe.sup[110]-0.015, "Support", cex=0.9, col=rgb(0,.4,.5,1))
text(x=110, y=pe.sat[110]+0.015, "Satisfaction", cex=0.9, col=rgb(1,.55,.3,1))
mtext(text="GDP growth increases by 1 SD in year 0", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Violent crime sim 

summary(supdem.91.t2$Viol_within)
sd(supdem.91.t2$Viol_within, na.rm=TRUE)
chg.viol = 0.2

# support simulation

pred.data.sup = data.frame(DemSup_m1=0, DemSup_m2=0, DemSat=0, 
                           Poly=0, Lib=0, GDP=0, 
                           GDPgr=0, Infl=0, HAQI=0, 
                           Viol=c(rep(0, n.burn), rep(chg.viol, n.yrs-n.burn)), Corr=0,
                           DemSat_m1=0, DemSat_m2=0, DemSup=0)

sup.mod.nam = c("DemSup_m1", "DemSup_m2", "DemSat", "Poly", "Lib", "GDP", "GDPgr", "Infl", "HAQI", "Viol", "Corr")
sat.mod.nam = c("DemSat_m1", "DemSat_m2", "DemSup", "Poly", "Lib", "GDP", "GDPgr", "Infl", "HAQI", "Viol", "Corr")

n.sims = 1e3
sim.data = array(as.matrix(pred.data.sup), dim=c(dim(pred.data.sup), n.sims))
mod.sims = mvtnorm::rmvnorm(n=n.sims, mean=coef(mod), sigma=vcov(mod))

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    sim.data[i, match("DemSat", names(pred.data.sup)), k] =
      mod.sims[k,(n.coef.sat+1):n.coef] %*% as.matrix(sim.data[i, match(sat.mod.nam, names(pred.data.sup)), k])
    sim.data[i+1, match("DemSat_m1", names(pred.data.sup)), k] = sim.data[i, match("DemSat", names(pred.data.sup)), k]
    sim.data[i+1, match("DemSat_m2", names(pred.data.sup)), k] = sim.data[i, match("DemSat_m1", names(pred.data.sup)), k]
    sim.data[i, match("DemSup", names(pred.data.sup)), k] = 
      mod.sims[k,1:n.coef.sup] %*% as.matrix(sim.data[i, match(sup.mod.nam, names(pred.data.sup)), k]) 
    sim.data[i+1, match("DemSup_m1", names(pred.data.sup)), k] = sim.data[i, match("DemSup", names(pred.data.sup)), k]
    sim.data[i+1, match("DemSup_m2", names(pred.data.sup)), k] = sim.data[i, match("DemSup_m1", names(pred.data.sup)), k]
  }
}

pe.sup = apply(sim.data[, match("DemSup", names(pred.data.sup)), ], 1, mean, na.rm=TRUE)
pe.sd = apply(sim.data[, match("DemSup", names(pred.data.sup)), ], 1, sd, na.rm=TRUE)
u95.sup = pe.sup + 1.96*pe.sd
l95.sup = pe.sup - 1.96*pe.sd

yr.st = n.burn - 1

# satisfaction simulation

pred.data.sat = data.frame(DemSup_m1=0, DemSup_m2=0, DemSat=0, 
                           Poly=0, Lib=0, GDP=0, 
                           GDPgr=0, Unemp=0, Infl=0, HAQI=0, 
                           Viol=c(rep(0, n.burn), rep(chg.viol, n.yrs-n.burn)), Corr=0,
                           DemSat_m1=0, DemSat_m2=0, DemSup=0)

n.sims = 1e3
sim.data = array(as.matrix(pred.data.sat), dim=c(dim(pred.data.sat), n.sims))

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    sim.data[i, match("DemSup", names(pred.data.sat)), k] = 
      mod.sims[k, 1:n.coef.sup] %*% as.matrix(sim.data[i, match(sup.mod.nam, names(pred.data.sat)), k]) 
    sim.data[i+1, match("DemSup_m1", names(pred.data.sat)), k] = sim.data[i, match("DemSup", names(pred.data.sat)), k]
    sim.data[i+1, match("DemSup_m2", names(pred.data.sat)), k] = sim.data[i, match("DemSup_m1", names(pred.data.sat)), k]
    sim.data[i, match("DemSat", names(pred.data.sat)), k] = 
      mod.sims[k,(n.coef.sat+1):n.coef] %*% as.matrix(sim.data[i, match(sat.mod.nam, names(pred.data.sat)), k]) 
    sim.data[i+1, match("DemSat_m1", names(pred.data.sat)), k] = sim.data[i, match("DemSat", names(pred.data.sat)), k]
    sim.data[i+1, match("DemSat_m2", names(pred.data.sat)), k] = sim.data[i, match("DemSat_m1", names(pred.data.sat)), k]
  }
}

pe.sat = apply(sim.data[, match("DemSat", names(pred.data.sat)), ], 1, mean, na.rm=TRUE)
pe.sd = apply(sim.data[, match("DemSat", names(pred.data.sat)), ], 1, sd, na.rm=TRUE)
u95.sat = pe.sat + 1.96*pe.sd
l95.sat = pe.sat - 1.96*pe.sd

# plot 

pdf("fig_3b.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe.sup[yr.st:(n.yrs-1)], type="n", ylim=c(-0.2, 0.2), 
     xaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(110, 120, by=10), lty=3, lwd=0.75, col=grey(0.5))
abline(h=c(-0.2, -0.1, 0, 0.1, 0.2), lty=3, lwd=0.75, col=grey(0.5))
abline(v=100, lty=2, lwd=0.75, col=grey(0.5))
axis(side=1, at=seq(100, 120, by=10), labels=seq(0, 20, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-0.2, -0.1, 0, 0.1, 0.2), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe.sup[yr.st:(n.yrs-1)], col=rgb(0,.4,.5,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(0,.4,.5,.5), border=NA, 
        y=c(u95.sup[yr.st:(n.yrs-1)], l95.sup[(n.yrs-1):yr.st]))
lines(x=yr.st:(n.yrs-1), y=pe.sat[yr.st:(n.yrs-1)], col=rgb(1,.55,0,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(1,.55,0,.5), border=NA, 
        y=c(u95.sat[yr.st:(n.yrs-1)], l95.sat[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Democratic support / satisfaction", cex=0.9)
text(x=110, y=pe.sup[110]+0.015, "Support", cex=0.9, col=rgb(0,.4,.5,1))
text(x=110, y=pe.sat[110]-0.015, "Satisfaction", cex=0.9, col=rgb(1,.55,.3,1))
mtext(text="Interpersonal violence increases by 1 SD in year 0", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Corruption effects sim 

mod = sem.fit.full.91
round(coef(sem.fit.full.91), 2)
n.coef.sup =  length(coef(mod)) / 2 - 1
n.coef.sat =  length(coef(mod)) / 2 - 1
n.coef = length(coef(mod)) - 2
vcov.mod = matrix(as.numeric(vcov(mod)), nrow=length(coef(mod)), ncol=length(coef(mod)))
n.yrs = 121
n.burn = 100

summary(supdem.91.t2$Corr_BCI_within)
sd(supdem.91.t2$Corr_BCI_within, na.rm=TRUE)
chg.corr = 0.14

# support simulation

pred.data.sup = data.frame(DemSup_m1=0, DemSup_m2=0, DemSat=0, 
                           Poly=0, Lib=0, GDP=0, 
                           GDPgr=0, Infl=0, HAQI=0, Viol=0, 
                           Corr=c(rep(0, n.burn), rep(chg.corr, n.yrs-n.burn)),
                           DemSat_m1=0, DemSat_m2=0, DemSup=0)

sup.mod.nam = c("DemSup_m1", "DemSup_m2", "DemSat", "Poly", "Lib", "GDP", "GDPgr", "Infl", "HAQI", "Viol", "Corr")
sat.mod.nam = c("DemSat_m1", "DemSat_m2", "DemSup", "Poly", "Lib", "GDP", "GDPgr", "Infl", "HAQI", "Viol", "Corr")

n.sims = 1e3
sim.data = array(as.matrix(pred.data.sup), dim=c(dim(pred.data.sup), n.sims))
mod.sims = mvtnorm::rmvnorm(n=n.sims, mean=coef(mod), sigma=vcov.mod)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    sim.data[i, match("DemSat", names(pred.data.sup)), k] =
      mod.sims[k,(n.coef.sat+1):n.coef] %*% as.matrix(sim.data[i, match(sat.mod.nam, names(pred.data.sup)), k])
    sim.data[i+1, match("DemSat_m1", names(pred.data.sup)), k] = sim.data[i, match("DemSat", names(pred.data.sup)), k]
    sim.data[i+1, match("DemSat_m2", names(pred.data.sup)), k] = sim.data[i, match("DemSat_m1", names(pred.data.sup)), k]
    sim.data[i, match("DemSup", names(pred.data.sup)), k] = 
      mod.sims[k,1:n.coef.sup] %*% as.matrix(sim.data[i, match(sup.mod.nam, names(pred.data.sup)), k]) 
    sim.data[i+1, match("DemSup_m1", names(pred.data.sup)), k] = sim.data[i, match("DemSup", names(pred.data.sup)), k]
    sim.data[i+1, match("DemSup_m2", names(pred.data.sup)), k] = sim.data[i, match("DemSup_m1", names(pred.data.sup)), k]
  }
}

pe.sup = apply(sim.data[, match("DemSup", names(pred.data.sup)), ], 1, mean, na.rm=TRUE)
pe.sd = apply(sim.data[, match("DemSup", names(pred.data.sup)), ], 1, sd, na.rm=TRUE)
u95.sup = pe.sup + 1.96*pe.sd
l95.sup = pe.sup - 1.96*pe.sd

yr.st = n.burn - 1

# satisfaction simulation

pred.data.sat = data.frame(DemSup_m1=0, DemSup_m2=0, DemSat=0, 
                           Poly=0, Lib=0, GDP=0, 
                           GDPgr=0, Infl=0, HAQI=0, Viol=0, 
                           Corr=c(rep(0, n.burn), rep(chg.corr, n.yrs-n.burn)),
                           DemSat_m1=0, DemSat_m2=0, DemSup=0)

n.sims = 1e3
sim.data = array(as.matrix(pred.data.sat), dim=c(dim(pred.data.sat), n.sims))

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    sim.data[i, match("DemSup", names(pred.data.sat)), k] = 
      mod.sims[k, 1:n.coef.sup] %*% as.matrix(sim.data[i, match(sup.mod.nam, names(pred.data.sat)), k]) 
    sim.data[i+1, match("DemSup_m1", names(pred.data.sat)), k] = sim.data[i, match("DemSup", names(pred.data.sat)), k]
    sim.data[i+1, match("DemSup_m2", names(pred.data.sat)), k] = sim.data[i, match("DemSup_m1", names(pred.data.sat)), k]
    sim.data[i, match("DemSat", names(pred.data.sat)), k] = 
      mod.sims[k,(n.coef.sat+1):n.coef] %*% as.matrix(sim.data[i, match(sat.mod.nam, names(pred.data.sat)), k]) 
    sim.data[i+1, match("DemSat_m1", names(pred.data.sat)), k] = sim.data[i, match("DemSat", names(pred.data.sat)), k]
    sim.data[i+1, match("DemSat_m2", names(pred.data.sat)), k] = sim.data[i, match("DemSat_m1", names(pred.data.sat)), k]
  }
}

pe.sat = apply(sim.data[, match("DemSat", names(pred.data.sat)), ], 1, mean, na.rm=TRUE)
pe.sd = apply(sim.data[, match("DemSat", names(pred.data.sat)), ], 1, sd, na.rm=TRUE)
u95.sat = pe.sat + 1.96*pe.sd
l95.sat = pe.sat - 1.96*pe.sd

# plot

pdf("fig_3c.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe.sup[yr.st:(n.yrs-1)], type="n", ylim=c(-0.2, 0.2), 
     xaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(110, 120, by=10), lty=3, lwd=0.75, col=grey(0.5))
abline(h=c(-0.2, -0.1, 0, 0.1, 0.2), lty=3, lwd=0.75, col=grey(0.5))
abline(v=100, lty=2, lwd=0.75, col=grey(0.5))
axis(side=1, at=seq(100, 120, by=10), labels=seq(0, 20, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-0.2, -0.1, 0, 0.1, 0.2), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe.sup[yr.st:(n.yrs-1)], col=rgb(0,.4,.5,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(0,.4,.5,.5), border=NA, 
        y=c(u95.sup[yr.st:(n.yrs-1)], l95.sup[(n.yrs-1):yr.st]))
lines(x=yr.st:(n.yrs-1), y=pe.sat[yr.st:(n.yrs-1)], col=rgb(1,.55,0,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(1,.55,0,.5), border=NA, 
        y=c(u95.sat[yr.st:(n.yrs-1)], l95.sat[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Democratic support / satisfaction", cex=0.9)
text(x=110, y=pe.sup[110]+0.015, "Support", cex=0.9,col=rgb(0,.4,.5,1))
text(x=110, y=pe.sat[110]+0.015, "Satisfaction", cex=0.9, col=rgb(1,.55,.3,1))
mtext(text="Corruption increases by 1 SD in year 0", side=3, line=0.4, cex=0.9, adj=0)
dev.off()

